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CO Automatic Decorrelation and Parameter Tuning Real-Time 
^ Kinematic Method and Apparatus 

Related Applications 

[001] The present application is related to co-pending U.S. Patent Application 
Serial Number 10/610,544 entitled, "Accurate and Fast RTK," and co-pending U.S. 
Patent Application Serial Number 10/610,541 entitled, "Enhanced Real Time 
Kinematics Determination Method and Apparatus," both assigned to the present 
assignee, both hereby incorporated herein by reference in their entirety, and both filed 
on July 2, 2003. 

Field of the Invention 

[002] The present invention relates to a method of and apparatus for an improved 
decorrelation and parameter tuning capabiHty for real time kinematics (RTK) 
determination. 

Background 

[003] High precision GPS kinematic positioning is widely used for many 
surveying and navigation applications on land, at sea and in the air. The distance 
from the mobile receiver to the nearest reference receiver may range from a few 
kilometers to hundreds of kilometers. As the receiver separation increases, the 
problem of accounting for distance-dependent biases grows and, as a consequence, 
reliable ambiguity resolution becomes an even greater challenge. 

[004] The standard mode of precise differential positioning uses one reference 
receiver located at a station whose coordinates are known, while determining a second 
receiver's coordinates relative to the reference receiver. In addition, the second 
receiver may be static or moving, and carrier phase measurements must be used to 
assure high positioning accuracy. This is the basis for pseudo-range-based differential 
global positioning system (GPS) (DGPS for short) techniques. However, for high 
precision applications, the use of carrier phase data comes at a cost in terms of overall 



system complexity because the measurements are ambiguous, requiring that 
ambiguity resolution (AR) algorithms be incorporated as an integral part of the data 
processing software. 

[005] Such high accuracy techniques result from progressive research and 
development (R&D) innovations, subsequently implemented by GPS manufacturers 
in top-of-the-line "GPS surveying" products. Over the last decade, several significant 
developments have resulted in the high accuracy performance also being available in 
"real-time" — that is, in the field, immediately following measurement, and after the 
data from the reference receiver has been received by the (second) field receiver for 
processing via a data communication link (e.g., very high frequency (VHP) or ultra 
high frequency (UHF) radio, cellular telephone, frequency modulation (FM) radio 
sub-carrier or satellite communication link). Real-time precise positioning is even 
possible when the GPS receiver is in motion through the use of "on-the-fly" (OTP) 
AR algorithms. These systems are commonly referred to as "real-time kinematic" 
(RTK) systems, and make feasible the use of GPS-RTK for many time-critical 
applications, e.g., machine control, GPS-guided .earthworks/excavations, automated 
haul truck operations, and other autonomous robotic navigation applications. 

[006J If the GPS signals were continuously tracked (loss-of-lock never occurred), 
the integer ambiguities resolved at the beginning of a survey would be valid for the 
whole GPS kinematic positioning span. However, the GPS sateUite signals are 
occasionally shaded (e.g. due to buildings in "urban canyon" environments), or 
momentarily blocked (e.g. when the receiver passes under a bridge or through a 
tunnel), and in these circumstances the integer ambiguity values are "lost" and must 
be re-determined or re-initialized. This process can take from a few tens of seconds up 
to several minutes with present OTF AR techniques. During this "re-initialization" 
period, the GPS carrier-range data cannot be obtained and there is "dead" time until 
sufficient data has been collected to resolve the ambiguities. If GPS signal 
interruptions occur repeatedly, ambiguity re-initialization is, at the very least, an 
irritation, and, at worst, a significant weakness of commercial GPS-RTK positioning 
systems. In addition, the longer the period of tracking required to ensure reliable OTP 
AR, the greater the risk that cycle slips occur during the crucial (re-)initialization 
period. A loss of lock of a receiver phase lock loop causing a sudden integer number 
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of cycles jump in a carrier phase observable is known as a cycle slip. Receiver 
tracking problems or an interrupted ability of the antenna to receive satellite signals 
causes the loss of lock. As the receiver separation increases, the problems of 
accounting for distance-dependent biases grow and, as a consequence, reliable 
ambiguity resolution (or re-initialization) becomes an even greater challenge due to: 

Residual biases or errors; 

Observation span determination; and 

Kinematic positioning results quality control. 

[007] Residual biases or errors after double-differencing can only be neglected 
for AR purposes when the distance between two receivers is less than about 10km. 
For distances greater than 10km, the distance-dependent biases, such as orbit bias, 
ionospheric delay and tropospheric delay, become significant problems. 

[008] Determining how long an observation span should be for reliable AR is a 
challenge for real-time GPS kinematic positioning. The longer the observation span 
required, the longer the "dead" time during which precise position determination is 
not possible. This can happen at the ambiguity initiaUzation step if a GPS survey is 
just starting, or at the ambiguity re-initialization step if the GPS signals are blocked, 
such as due to cycle slips or data interruptions. 

[009] Quality control of GPS kinematic positioning results is a critical issue and 
is necessary during all processes: data collection, data processing and data 
transmission. Quality control procedures are not only applied for carrier phase-based 
GPS kinematic positioning, but also for pseudo-range-based DGPS positioning. 
However, quality control or validation criteria for AR, for precise GPS kinematic 
positioning, is a significant challenge. 

[010] There is a need in the art for an improved method of processing GPS 
signals to improve decorrelation and parameter tuning capabilities for high accuracy, 
reUable position determination. 
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Summary 

[Oil] It is therefore an object of the present invention to provide a method of 
processing GPS signals to improve decorrelation and parameter tuning capabilities for 
high accuracy, reliable position determination. 

[012] In an embodiment according to the present application, the ADAPT-RTK 
Technology (Automatic Decorrelation and Parameter Tuning Real-Time Kinematic 
technology) is described using an innovative decorrelation algorithm to decorrelate 
the correlation between LI and L2 measurements through a functional model and a 
stochastic model reducing the impact on performance of distance dependent biases. 
Therefore, the ADAPT-RTK Technology improves GPS RTK performance to exceed 
other products in the market, especially on long-range applications. 

[013] The present inventors have successfully implemented the present 
decorrelation approach in a Z-Max RTK engine from Thales Navigation. After 
extensive testing using Z-Xtreme and Z-Max receivers, the most exciting results have 
been derived. 

[014] The decorrelation approach has significant advantages: (1) it more closely 
describes the measurement natural features, especially in long-range applications; (2) 
it is a rigorous approach from a Kalman filtering implementation point of view; (3) all 
other approaches, e.g., widelane/iono-free approach, widelane only approach, become 
special cases in theory and in practice, and; (4) it significantly improves RTK 
performance, especially for long-range applications. 

[015] The functional model includes determining the reference satellite, residuals, 
and design matrix. The stochastic model includes variance-covariance matrix 
determination for measurements and the variance-covariance matrix for dynamic 
noise in Kalman filtering. 

[016] The RTK technique requires available dual frequency pseudo-range and 
carrier phase measurements. On-the-fly RTK might not necessarily need pseudo- 
range measurements because the float solution can be derived by the change of carrier 
phase measurements between epochs (or Doppler measurements). However, the float 
solution using data from a single epoch must be derived by pseudo-range 
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measurements. The integrated function model means modeling carrier phase 
measurements and also pseudo-range measurements and their stochastic features are 
included. The stochastic model for carrier phase measurements and pseudo-range 
measurements, especially the ratio between standard deviations of carrier phase and 
pseudo-range, significantly affects the float solution and subsequently affects the 
RTK performance. The present invention provides appropriate stochastic models for 
both carrier phase and pseudo-range to enhance RTK performance. 

[017] Still other objects and advantages of the present invention will become 
readily apparent to those skilled in the art from the following detailed description, 
wherein the preferred embodiments of the invention are shown and described, simply 
by way of illustration of the best mode contemplated of carrying out the invention. 
As will be realized, the invention is capable of other and different embodiments, and 
its several details are capable of modifications in various obvious respects, all without 
departing from the invention. 

Description of the Drawings 

[018] The present invention is illustrated by way of example, and not by 
limitation, in the figures of the accompanying drawings, wherein elements having the 
same reference numeral designations represent like elements throughout and wherein: 

Figure 1 is a graph of coefficients for a correlation approach according to an 
embodiment of the present invention; 

Figure 2 is a flow diagram of a decorrelation approach according to an 
embodiment of the present invention; 

Figure 3 is a graph of standard deviation for carrier phase and pseudorange 
according to an embodiment of the present invention; 

Figure 4 is a high level flow diagram of a process according to an embodiment 
of the present invention; 

Figure 5 is a flow diagram of a portion of the Figure 4 flow diagram according 
to an embodiment of the present invention; and 
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Figure 6 is a computer system on which an embodiment of the present 
invention may be used. 

Detailed Description 

[019] The functional model of the present invention includes a Double 
Differenced Model. The double-differenced observable in units of meters can be 
formed as: 

if the four reference satellites for GPS LI, GPS L2, GLONASS LI and GLONASS 
L2 are chosen. In this case, four receiver clock term and inter-channel biases can be 
removed. 

[020] The integer ambiguity terms can also be re-arranged in the following way: 
^ j^* j - A<|)k = VApk^j +XyVAU + (X J- - ). AN,, + Al/f ? - Al/f^ + 8 va*,, ( 1 -2) 

For GPS carrier phase measurements, the third term on the lefthand side of equation 
(1-2) disappears because Aj and are equal. For GLONASS carrier phase 
measurements, the third term, or the single-differenced integer ambiguity for the 
reference satellite, must be estimated prior to estimating the double-differenced 
integer ambiguities. The remaining errors from the third term can cause systematic 
model errors and may result in incorrect double-differenced ambiguity resolution, and 
hence degraded positioning accuracy. 

[021] An alternative approach is to form the double-differenced observable after 
the single-differenced observables are expressed in units of cycles as in equation 1-3, 
as follows: 




(1-3) 

where AdT^ are different for GPS LI, GPS L2, GLONASS LI and GLONASS L2 
signals. The differenced receiver clock bias (or inter-channel biases) cannot be 
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eliminated in equation (1-3). Tlie fourth term (ionospheric delay) and the fifth term 
(tropospheric delay) become slightly larger than in the case when the two frequencies 
are the same. Using GPS and GLONASS pseudo-range measurements, the difference 
between the two receiver clock biases is estimated. The estimate could be used to 
correct the second term for ambiguity resolution purposes; however, the receiver 
clock bias significantly degrades positioning accuracy. 

[022] In an embodiment according to the present invention, the proposed 
approach estimates AdT^ using pseudorange and the second term is corrected prior to 
ambiguity resolution by using equation (1-3). After fixing the ambiguity, equation (1- 
2) is used for position fixing where the third term on the right side of equation (1-2) is 
estimated. For a GPS-only positioning system, the. third term on the right side of 
equation (1-2) is always zero because the frequencies for satellite j and k are the same. 
The ionosphere delay is corrected by a GPS broadcast model (Klobuchar Model). 
The troposphere delay is corrected by the Hopfield model known to persons of skill in 
the art. 

Kalman Filtering State Vector 

[023] For a GPS-only RTK system, the following states in a Kalman filter have 
been implemented: 



Elements in State Vector 


Dimension 


Notes 


Position component X (or Easting) 




Mandatory 


Velocity component X (or Easting) 




Optional 


Acceleration component X (or Easting) 




Optional 


Position component Y (or Northing) 




Mandatory 


Velocity component Y (or Northing) 




Optional 


Acceleration component Y (or Northing) 




Optional 


Position component Z (or Height) 




Mandatory 


Velocity component Z (or Height) 




Optional 


Acceleration component Z (or Height) 




Optional 


Zenith troposphere delay scale factor 




Optional 


Residual DD ionosphere delay 


Nsat-1 


Optional 


LI DD ambiguity 


Nsat-1 


Optional 


L2 DD ambiguity 


Nsat-1 


Optional 



Notes: Nsat is the number of visible satellites. 
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[024] The ionosphere delay and the zenith troposphere delay scale factor are 
estimated after model correction although they operate without model corrections. 

Transition Matrix and Dynamic Noise 

Position Parameters 

[025] When an observer is nearly stationary, such as a buoy drifting at sea, or 
crustal deformation monitoring, the position is assumed to be a random-walk process. 
In this case, three coordinate parameters are enough to employ the Kalman state 
vector. The transition matrix and dynamic noise are determined based on random 
walk model, such as equations 1-4 and 1-5: 



«t>M-.=e'^*'""-'^=l (1-4) 
Q. =<^u(t. -t,..) (1-5) 



[026] 



which is called the Position model. 



[027] When the observer is not stationary, but moving with nearly constant 
velocity, the velocity is not white noise but a random-walk process. In this case, three 
coordinate parameters and three velocity parameters must be included in the Kalman 
state vector. The transition matrix and dynamic noise can be determined based on 
integrated random walk model, such as equations 1-6 and 1-7: 



[028] 



<l>k,k-i = 



1 ih-hJ 
0 1 



(1-6) 



[029] 



Qk=^u 



(1-7) 



[030] 



which is called the Position-Velocity model. 



[031] The Position-Velocity model also becomes inadequate for cases where the 
near-constant velocity assumption is incorrect, that is, in the presence of severe 
accelerations. To accommodate acceleration in the process model, it is appropriate to 
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add another degree of freedom for each position state which becomes Position- 
Velocity-Acceleration model or Gauss-Markov process is perhaps more appropriate 
than the non-stationary random walk for acceleration. 

[032] Determination of the spectral amplitude for position random processes is a 
best "guess" based on expected vehicle dynamics. In many vehicular applications, 
random perturbations are greater in the horizontal plane than in the vertical plane. 
These perturbations are accounted for by selecting a value for spectral amplitude 
which is lesser for the altitude channel, i.e., vertical plane, than for the two horizontal 
ones, i.e., horizontal plane. 

Troposphere Delay 

[033] The troposphere delay includes a wet component and a dry component. 
Whereas the dry component can be modeled with a high degree of accuracy, the wet 
component cannot be so modeled. A solution is to scale the wet component by a scale 
factor, i.e., (l + e) . After troposphere delay modeling correction, the estimation part 
is e-Ad^^^^f. Therefore, the design matrix element is Adf^^ and e represents the 

changing percentage of the wet component from the modeled value. Satellite 
elevation has no impact as the wet component for all satellites is scaled by the same 
factor for a given location. Empirically, the scale factor is modeled as a Gauss- 
Markov process. The transition matrix and dynamic model is then derived as in 
equations 1-8 and 1-9, as follows: 

*k.k-i ^e"^*-^*^"*^-'^ (1-8) 

Qk=a^^p(l-e-'P-^^*^-*-'>) (1-9) 

[034] where l/Pt^^p the correlation time of the troposphere wet component and 
^Lp represents the wet component changing level, both of which are functions of the 

baseline length and the height difference. The correlation time is set based on 
empirical data. The default value is set at 600 seconds. The initial variance for zenith 
troposphere scale factor is set as a default value of 0.1 (unitless) and dynamic noise is 
set as a default value of 0.1 ppm. 
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Ionosphere delay Parameters 

[035] The design matrix element of an ionosphere delay parameter is a mapping 
function of the satellite elevation and the difference of the ionosphere zenith delays on 
two intersections between the receiver-satellite rays from both ends of a baseline and 
an equivalent ionosphere layer. However, a part of the ionosphere delay is also 
corrected based on the ionosphere delay model. 

[036] The ionosphere zenith delay, Ad?„, is a function of the local time, 
ionosphere activities, distance and direction of two intersections of receiver-satellite 
rays with equivalent ionosphere layer from both ends of a baseline. The parameters 
of the empirical model to approximate an ionosphere delay have been broadcast by 
GPS satellites in real-time, which can be used for compensating for ionosphere delay. 
This model provides at least a 50% reduction in the single-frequency user's root mean 
square (RMS) due to ionosphere propagation effects. The model is widely used in a 
high-precision carrier phase-based RTK engine. The remaining part of the 
ionosphere delay is estimated in Kalman filtering as a state vector element and is 
empirically estimated using a Gauss-Markov model. The transition matrix and 
dynamic model are expressed using equations 1-10 and 1-1 1, as follows: 

<t>k,k-i =e-''^(*^-''=-'' (1-10) 

Qk=CTLO-e-^P-^*^-*-^) (1-11) 

[037] where l/P-^^ is the correlation time of Adfo„ after ionosphere broadcast 
model corrections, and a?n represents variation level of the Ad? „ . 

[038] The initial standard deviation for residual ionosphere parameters is set to a 
default value of 0.5ppm*distance. The dynamic noise for residual ionosphere 
parameters is 0.01ppm*sqrt(distance) and the weighting function for the different 
satellite elevation angle is 1 .0/sqrt(sqrt(cos(90-elev))). 

Decorrelation Approach and Implementation 

Decorrelation 
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[046] The measurement and modeling errors consist of measurement noise, 
multipath, ionosphere delay, troposphere delay, orbit bias, inter-channel bias and 
antenna offset and additional error sources indicated from warning messages. These 
errors, or biases, are classified into two categories, e.g. distance-dependent biases and 
distance-independent biases, and represented by equations 2-1 and 2-2 as follows: 

^non-dist = K-Lise + ^MP + ^Im (2-1) 

RL=RL+RLp+Rorb (2-2) 

[047] The errors from inter-channel bias, antenna offset etc. are assumed to be 
cancelled through double difference. Only the errors remaining in double differenced 
measurements are taken into account. 

[048] Due to the distance-dependent errors such as ionosphere, troposphere and 
orbit bias, the LI and L2 carrier phase measurements become highly correlated when 
baseline length increases and the errors are not processed independently as performed 
in the short range case. Therefore, the co-variance matrix between LI and L2 must be 
taken into account. The detailed implementation is as follows in equations 2-3, 2-4, 
and 2-5: 

[049] V, =H,X + L, (2-3) 
[050] V3=H,X + L2 (2-4) 

(2-5) 

[052] where is pre-fit residuals for the LI measurement vector (the difference 
between the computed vector and the measurement vector), is a design matrix, ^' 
is a variance-covariance matrix of the LI measurement vector, and X represents the 
Kalman filter state vector, is a post-fit residual vector. ^2, ^2, ^2 have 

similar meanings for the L2 measurement vector. ^12 or ^21 is a co-variance matrix 
between the LI and L2 measurement vectors. 



[051] 



D(L„L,) = 
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[046] Due to the correlation of the LI and L2 measurement vectors, LI block and 
L2 block are not updated independently. The correlation approach orthogonalizes the 

L2 measurement vector and forms a new measurement vector ^ which is 
independent from the LI measurement vector. The new measurement is as follows in 
equation 2-6: 



[047] where: 



=H2X + L2 



(2-6) 



Lj =L2 -D^.D^'L, (2-7) 



H,=H,-D„Dr'H, (2-8) 



-Dj,D-'D,2, and (2-9) 



D(L„L,) = 



D, 
0 



0 



(2-10) 



T FT 
[048] After the new measurement vector 2^ design matrix 2 and variance 

matrix ^2 are formed, LI block and ^2 block are updated independently. However, 
there is some additional computation load, as described in detail below. 



Implementation 

R ^ o 2 

[049] The zenith variances li^ ^L2for single differenced (between receivers) 
LI and L2 measurements consist of non-distance dependent errors and distance 
dependent errors, e.g., ionosphere delay, troposphere delay, and orbit bias, as follows 
in equations 2-11 and 2-12: 



^L2 -TT 



R ^ 4- 



V f 



ion ^trop ^orb 



(2-11) 



(2-12) 
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[058] where B is the baseline length. The correlated errors " "between LI and 
L2 zenith is as follows in equation 2-13: 



(2-13) 



[059] The variance at other elevation angles is mapped by a function, e.g., 
equation 2-14: 



W = 



fl.0 + 2.5e'^"^ forpseudorange 
[l.0 + 7.5-e"^"^ for carrier phase 



(2-14) 



[060] The variance and co-variance matrix for the single difference measurement 
vectors is represented as follows using equation 2-15: ' 



where 



SD = 



SD,, SD,, 
SD^i SD^ij (2-15) 



SD22 = 2K\ 



w, 


0 • 


0 


0 ■ 


0 


w, • 


0 


0 


0 


0 • 




0 


0 


0 


0 


w„. 


w, 


0 


0 


0 ' 


0 




0 


0 


0 


0 


•• w„_, 


0 


0 


0 


• 0 


Wo. 



(2-16) 



(2-17) 
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SD,2 - SDj, - 2Rl,_l2 



W, 0 

0 

0 0 

0 0 



0 0 

0 0 

W„.. 0 

0 W„ 



(2-18) 



[061] The variance and co-variance matrix for double difference measurements 
(between receivers and satellites) is represented as follows using equation 2-19: 



DD = 



(2-19) 



where 



DD„ =2R^, 



w„ 



w. 



w„ 



w. 



ref 



w, 



ref 



w. 



ref 



w. 



ref 



w + w 

^ ^^ref 



W„ 



W„_, +W, 



ref. 



(2-20) 



DE>«=2RL 



w, + w„, w, 



w„ 



w„ 



ref ^'^ref 

w. 



w. 



ref 



w. 



ref 



w„ 

w, 



ref 



w. 



ref 



w„_,+w„ 



(2-21) 



DD„ = DD,3 = 2R„,,, 



'ref ^^ref 
^ref W, 



w, 
w, 



ref 



w„ 



ref 



w, 



ref 



w, 
w. 



ref 



ref 



ref 



(2-22) 



[062] Equation 2-23 for calculating ^^n^^n is; 
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dDijDd;^,' = 





"1 


0 • 


• 0 


0 


2Rl1,L2 


0 


1 


• 0 


0 




0 


0 • 


• 1 


0 




0 


0 • 


• 0 


1 



(2-23) 

[055] Hence, computation is simplified to a coefficient associated with a preset 
setting and baseline length as in equation 2-24 as follows: 



Coef=^^'-'" 



(2 2 2 /f,' 2 2 2 1 

(Rfon + Rnop + Roib J ^Rfon + ^Lp + Rprb 



^non-dist,LI ' 



> 2 
"•ion 



2 



(2-24) 



[056] There are three special cases for the coefficient factor, each of which is 
described in detail below. 

Case 1 : Ionosphere-Free combination 

[057] Assuming there is no non-distance dependent error, no troposphere delay 

13 2 T^2 p 2 

and no orbit bias (^non-dist.Li=o, »^°p=0, °^b=o), Coef of equation 2-24 is as 
follows: 



R 



Coef =- 



X.4 

^2 



-J-=-^ = 1.283 



A*2 ^1 



[058] This means the orthogonalized measurement (L2(cycle)- ^ *L1 (cycle)) is 
the ionosphere free combination. 

Case 2: LI and L2 Independent Measurements 

R^ 

[059] Assuming the distance-dependent biases are ignored, ( =0, =0 and 



R 



*™p=0), Coef of equation 2-24 is as follows: 



Coef - 0 
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[060] This means the orthogonalized measurement is the L2 measurement. 
Case 3: Geometry-Free combination 

[061] Assuming the non-distance dependent error and ionosphere delay are zero 

( non-dist.Li =0, non-dist,L2 ^Q, and =0), frequcncy independent errors from 
distance-dependent errors such orbit errors and troposphere delay make up the 
dominator. Coef of equation 2-24 is as follows: 

Coef =-^ = 0.7792 

[062] This means the orthogonalized measurement will be L2(cycle)- ^ 
*L1 (cycle), which is the geometry free combination. 

[063] The decorrelation approach has been implemented in a Z-Xtreme-type and 
Z-Max-type receiver and the inventors have observed that the decorrelation approach 
selects the best LI and L2 combination according to preset settings. Based on the 
decorrelation approach, any combination between LI and L2 is easily implemented as 
follows from equation 2-25: 

[064] The new combination is formed as: 

¥2=^2^ -^^2 (2-25) 

where 

^2 =^2 (2-26) 
H^^H^-aH^^ (2-27) 
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D, =(R^,+a^R^„-2aRu.L2) 



w, 



ref 



w, 
w, 



ref 



ref 



w, 
w. 



ref 



ref 



(2-28) 



W. 



tef 



w. 



ref 



w, 
w, 



ref 



ref 



w„., + w^ 



[065] However, the correlation between LI and ^2 has not been taken into 
account although all combination options are implemented based on the subroutine 
for correlation approach. 



a = < 



1 

^2 Ai 
0 

Rl1,L2/R 



widelane combination 
geometry - free combination 
ionosphere - free combination 
L2 independent measurement 
the correlation approach 



(2-29) 



[066] The coefficients for the decorrelation approach have been plotted in Figure 

1 . The Kalman filtering update for the decorrelation approach is illustrated in Figure 

2, wherein the switch to determine ^ is determined based on equation (2-29). 

Stochastic Model 

[067] The measurement and modeling errors consist of measurement noise, 
multipath, ionosphere delay, troposphere delay and orbit bias as described by 
equations 2-11 and 2-12 above. 

[068] The remaining errors based on the ionosphere delay, troposphere and orbit 
bias are represented by distance-dependent functions for their standard deviations, eg. 
O.Sppm, 10-4 * HeightDiff + 2ppm, and O.lppm, respectively. If the ionosphere 
parameters are estimated in the Kalman filtering state vector, the standard deviation 
for ionosphere delay is scaled by a factor, e.g., 1 - dt/dTIono, where dt is the time 
shift from filtering start, dTIono is the time to turn off ionosphere. If the troposphere 
scale parameter is estimated, or a default standard troposphere model is applied, 
troposphere errors are turned off. 
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[069] By default, ionosphere (one for each Sat) and troposphere (one scale 
parameter) are modeled in Kalman filtering. In correspondence with the above 
paragraph, the troposphere term is not applied to the stochastic model. The 
ionosphere term is applied to the stochastic model even when the ionosphere is 
modeled. There is not any theoretical foundation behind the formula 1-dt/dTiono; this 
solution is heuristic and the inventors have found that the formula provides good 
results. Here, dt is the elapsed time fi-om starting estimation given ionosphere, dTiono 
is set between from 200 seconds to 1,000 seconds. The default value is set as 600 
seconds. 

[070] The standard deviation for constant errors from multipath and noise 
depends on a multipath template index. The index is defined as a value from 1 to 3 
representing the multipath status (1= minor multipath, 2 = medium, and 3 = strong 
multipath). 

[071] The default standard deviation for medium multipath have a default value 
of 2.25 meters for pseudo-range and 0.05 cycles for carrier phase. The mapping 
function is determined by the empirical exponential fiinction of equation 2-14. 

[072] The standard deviation for non-distance dependent errors is depicted in 
Figure 3 when the zenith default standard deviations are set at 0.05 cycle for phase 
and 2.25 meters for code. 

Ambiguity Resolution Validation 

[073] The ratio smoothing algorithm and ambiguity stability function are only 
used for the long range widelane/ionosphere-fi'ee combination. Both fimctions 
increase positioning reliability and increase time-to-first-fix. 

Instant Ambiguity Fixing 

[074] The methods claimed in the pending patent applications entitled, "Accurate 
and Fast RTK" and "Enhanced Real Time Kinematics Determination Method and 
Apparatus" and having application serial numbers 10/610,544 and 10/610,541, 
respectively, have been used for instant ambiguity fixing for baselines shorter than 10 
km. 
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Ratio Smoothing Algorithm for Long Range Applications 

[075] Ratio = Ration Ratio *(N--1)/N ^^^^^ ^j^^ ^j^^^ ^j^^ ^^^.^^ 

window (1 second per km, if the baseline length is longer than 10km). Both Ratio 
and Ratio ^ described in more detail in co-pending U.S. Patent Application Serial 
Numbers 10/610,544 and 10/610,541, are required to pass validation tests for 
ambiguity fixing. 

[076] A ratio smoothing algorithm increases reliability for long range; however, 
the same algorithm decreases instant-fix significantly. 

Ambiguity Repeatability Function for long range 

[077] The ambiguity repeatability check is for determining whether all double- 
difference ambiguities remain the same for N continuous epochs, e.g., 30, and 
whether all passed the ratio test mentioned above. If the ambiguity repeatability 
check is passed, the ambiguity is fixed. This function has a problem when cycle slips 
occur. The ambiguity stability counter is defined as a counter to count how many 
continuous epochs with the same ambiguity set or another ambiguity set with a 
difference larger than 50 cycles. This means that the ambiguity stability counter is 
not reset if the ambiguity difference is larger than 50 cycles. In this manner, the 
influence of cycle slips is minimized. This is reasonable because double differenced 
(DD) ambiguity is impossibly fixed to a wrong integer with such large difference 
values. However, the DD ambiguity stability counter is reset to a value of 1 if the 
ambiguity difference is not zero and less than 50 cycles for any given satellite pair. 

[078] The Ambiguity Repeatability Check increases reliability for long range 
applications; however, the same function makes instant-fix impossible. 

Partial Fixing Scheme 

[079] The partial fixing scheme has been implemented based on the functional 
model improvement by removing measurements and ambiguity searching adaptation 
by removing ambiguity, as described in co-pending U.S. Patent Application Serial 
Numbers 10/610,544 and 10/610,541. 
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Process flow 



[080] Figure 4 depicts the process flow of an RTK method according to an 
embodiment of the present invention. Low rate, nominally 1 Hz, base station 
measurements output from a base data decoder 400 are output to a polynomial fitting 
function executed in a phase predictor process 402, e.g. second order or higher order 
polynomial, and a Kalman filter process 408. Base data decoder 400 decodes raw 
GPS measurements received from a base GPS receiver (not shown). The sampling 
rate (or update rate) of base data decoder 400 is nominally 1 Hz. Higher sampling 
rates may be used with a proportional increase in cost and baud rate for the data link. 
In one embodiment of the present invention, the update rate does not exceed 1 Hz. If 
a higher sampling rate is required, an embodiment according to the description 
embodied in co-pending patent application titled, "Enhanced Rapid Real-Time 
Kinematic Determination Method and Apparatus," by the present inventors and 
assigned to the present assignee would be used. 

[081] In order to reduce the position update time delay, phase predictor process 
402 predicts corrections for the position calculation using available corrections 
transmitted from the base GPS receiver as decoded and output from base data decoder 
400 in conjunction with polynomial filtering. The positioning accuracy degrades 
depending on the length of the predicted period. 

[082] Kalman filter process 408, described in detail below with reference to 
Figure 5, calculates optimal solutions, i.e. position and/or velocity, based on currently 
available measurements from base data decoder 400 and rover data decoder 404. An 
ambiguity resolution process 410 is a part of Kalman filter process 408 and is 
described in detail below with reference to Figure 5. 

[083] A rover data decoder 404 decodes raw GPS measurements and ephemeris 
received from the rover GPS receiver (not shown) and provides time tagged carrier 
phase measurements to a carrier phase process 406. The sampling rate (or update 
rate) of rover data decoder 404 can be up to 10 Hz or higher. 
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[0102] Output estimated polynomial parameters from phase predictor 402 are used 
by carrier phase process 406 to predict the base station measurement at a rate 
matching the RTK update rate, typically 10 Hz or higher. The RTK solution latency 
is primarily determined by the rover measurement data collection time and the RTK 
position computation time. The base station prediction time is a negligible delay. 
Typically, RTK solution latency is less than 20 milliseconds depending on 
microprocessor speed. 

[0103] For embodiments requiring an update rate in time-tagged mode of 1 Hz or 
lower, carrier phase process 406 uses the output from Kalman filter 408 to calculate 
and output the rover GPS receiver position and/or velocity. For embodiments 
requiring an update rate in RTK mode of 1 Hz or lower, carrier phase process 406 
uses the most recent measurements from rover data decoder 404 and the predicted 
correction output from phase predictor 402 to calculate and output the most recent 
position and/or velocity of the rover GPS receiver. 

[0104] To reduce the RTK position computation time, only an LI carrier phase 
measurement output from a rover data decoder 404 is used. The rover LI carrier 
phase measurement and the predicted base station LI carrier phase measurement 
output from phase predictor 402 are then used to derive LI double difference 
measurements in 406. The estimated LI integer ambiguities, residual ionosphere 
delay, residual troposphere delay and other bias parameters are used to correct the 
double difference measurement in 406. The corrected double difference measurement 
is output into a least squares (LSQ) estimator to calculate a rover position in 406. The 
velocity is calculated in a similar manner using rover LI Doppler measurements and 
predicted base LI carrier phase rate. Because of the requirement of base station 
measurement prediction, the RTK solution accuracy is degraded in comparison to a 
matched time-tag RTK solution. With Selective Availability (S/A), the rate of 
degradation increases due to the inability to predict S/A. Selective availability is 
known to persons of skill in the art and refers to the intentional degradation of the 
absolute positioning performance capabilities of the GPS for civilian use 
accomplished by artificial "dithering" of satellite clock error. 

[0105] With reference to Figure 5, details of Kalman filter process 408 and 
ambiguity resolution process 410 in Figure 4 are now described. Data from both a 
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base GPS receiver (not shown) and a rover GPS receiver (not shown) is received, 
decoded, and output by base data decoder 400 and rover data decoder 404, 
respectively. Base data 500 and rover data 502 time tags are matched in match time 
tag step 504, thereby matching the time when the respective measurements were 
made. 

[0106] After the time tags are matched, the matched data output from match time 
tag step 504 is input to a Kalman filter. In step 506, if the matched data output is in 
the first epoch or if a reset of the Kalman filter is required, the Kalman filter is 
initialized in step 506. In step 508, a reference satellite is selected to determine the 
double differenced measurement. Further, cycle slips are checked using cycle slip 
flags and a stochastic model is calculated. 

[01071 The flow proceeds to step 510 for tiie preparation of the design matrix, 
variance-covariance matrix (stochastic model) and calculation of pre-fit residuals for 
all measurements, e.g. C/A pseudo-range, PI pseudo-range, P2 pseudo-range, LI 
Doppler, L2 Doppler, LI carrier phase and L2 carrier phase measurements. The 
output of the pre-fit residual calculation is input to a Receiver Autonomous Integrity 
Monitoring (RAIM) algorithm to detect outHers. RAIM is a form of receiver self- 
checking using redundant pseudo-range observations to detect if a problem with any 
of the measurements exists. 

[0108] The output of step 510 is provided as input to a Kalman filter measurement 
update step 512 to sequentially filter all measurements and provide filtered output 
measurements to an ambiguity resolution step 514. Figure 2 illustrates the detailed 
procedure implementing the decorrelation approach in Step 512. The update step 
provides the optimal estimation results using all available measurements. 

[0109] The validation criteria are calculated using the above-described method in 
step 514 and a determination of whether the integer ambiguities can be fixed or not is 
performed in step 516. If the step 516 determination is positive (the integer 
ambiguities can be fixed), the float solution is updated to the fix solution in step 518. 
If the step 516 determination is negative (the integer ambiguities cannot be fixed), the 
above-described adaptive fix procedure is performed in step 520 to attempt to fix 
ambiguities and a second determination of whether the integer ambiguities can be 
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fixed or not is performed in step 522. If the step 522 determination is positive (the 
integer ambiguities can be fixed), the float solution is updated to the fix solution in 
step 518 and the flow proceeds to step 524 wherein the post-fit residuals are updated 
and possible outliers are detected. If the step 522 determination is negative (the 
integer ambiguities cannot be fixed), the flow proceeds to step 524 described above. 

[0110] The output of step 524 is provided to a Kalman filtering time update 526, 
which is a Kalman filtering prediction step. In step 528, all necessary information is 
outputted and in step 530 the measurements are stored and the processing information 
is updated based on the above-described method. The flow proceeds to process the 
next epoch of data returning to step 504. 

[Dili] In coordination with the above-described technique, an embodiment of the 
present invention provides an improved method of and apparatus for an improved 
decorrelation and parameter tuning capability for determining real time kinematics. 

[0112] Figure 6 is a block diagram illustrating an exemplary computer 600 upon 
which an embodiment of the invention may be implemented. The present invention is 
usable with currently available handheld and embedded devices, e.g. GPS receivers, 
and is also applicable to personal computers, mini-mainfi:ames, servers and the like. 

[0113] Computer 600 includes a bus 602 or other communication mechanism for 
communicating information, and a processor 604 coupled with the bus 602 for 
processing information. Computer 600 also includes a main memory 606, such as a 
random access memory (RAM) or other dynamic storage device, coupled to the bus 
602 for storing GPS data signals according to an embodiment of the present invention 
and instructions to be executed by processor 604. Main memory 606 also may be 
used for storing temporary variables or other intermediate information during 
execution of instructions to be executed by processor 604. Computer 600 fiirther 
includes a read only memory (ROM) 608 or other static storage device coupled to the 
bus 602 for storing static information and instructions for the processor 604. A 
storage device 610 (dotted line), such as a compact flash, smart media, or other 
storage device, is optionally provided and coupled to the bus 602 for storing 
instructions. 
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[0114] Computer system 600 may be coupled via the bus 602 to a display 612, 
such as a cathode ray tube (CRT) or a flat panel display, for displaying an interface to 
the user. An input device 614, including alphanumeric and function keys, is coupled 
to the bus 602 for communicating information and command selections to the 
processor 604. Another type of user input device is cursor control 616, such as a 
mouse, a trackball, or cursor direction keys for communicating direction information 
and command selections to processor 604 and for controlling cursor movement on the 
display 612. This input device typically has two degrees of freedom in two axes, a 
first axes (e.g., x) and a second axis (e.g., y) allowing the device to specify positions 
in a plane. 

[0115] The invention is related to the use of computer 600, such as the depicted 
computer of Figure 6, to perform accurate, real-time, kinematics determination. 
According to one embodiment of the invention, data signals are received via a 
navigation interface 619, e.g. a GPS receiver, and processed by computer 600 aiid 
processor 604 executes sequences of instructions contained in main memory 606 in 
response to input received via input device 614, cursor control 616, or communication 
interface 618. Such instructions may be read into main memory 606 from another 
computer-readable medium, such as storage device 610. A user interacts with the 
system via an application providing a user interface displayed on display 612. 

[0116] However, the computer-readable medium is not limited to devices such as 
storage device 610. For example, the computer-readable medium may include a 
floppy disk, a flexible disk, hard disk, magnetic tape, or any other magnetic medium, 
a compact disc-read only memory (CD-ROM), any other optical medium, punch 
cards, paper tape, any other physical medium with patterns of holes, a random access 
memory (RAM), a programmable read only memory (PROM), an erasable PROM 
(EPROM), a Flash-EPROM, any other memory chip or cartridge, a carrier wave 
embodied in an electrical, electromagnetic, infrared, or optical signal, or any other 
medium from which a computer can read. Execution of the sequences of instructions 
contained in the main memory 606 causes the processor 604 to perform the process 
steps described above. In alternative embodiments, hard-wired circuitry may be used 
in place of or in combination with computer software instructions to implement the 
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invention. Thus, embodiments of the invention are not limited to any specific 
combination of hardware circuitry and software. 

[0117] Computer 600 also includes a communication interface 618 coupled to the 
bus 602 and providing two-way data communication as is known in the art. For 
example, communication interface 618 may be an integrated services digital network 
(ISDN) card, a digital subscriber line (DSL) card, or a modem to provide a data 
communication connection to a corresponding type of telephone line. As another 
example, communication interface 618 may be a local area network (LAN) card to 
provide a data communication coimection to a compatible LAN. Wireless links may 
also be implemented. In any such implementation, communication interface 618 
sends and receives electrical, electromagnetic or optical signals which carry digital 
data streams representing various types of information. Of particular note, the 
communications through interface 618 may permit transmission or receipt of 
instructions and data to be processed according to the above method. For example, 
two or more computers 600 may be networked together in a conventional manner 
with each using the communication interface 618. 

[0118] Network link 620 typically provides data communication through one or 
more networks to other data devices. For example, network link 620 may provide a 
connection through local network 622 to a host computer 624 or to data equipment 
operated by an Internet Service Provider (ISP) 626. ISP 626 in turn provides data 
communication services through the world wide packet data communication network 
now commonly referred to as the "Internet" 628. Local network 622 and Internet 628 
both use electrical, electromagnetic or optical signals which carry digital data streams. 
The signals through the various networks and the signals on network link 620 and 
through communication interface 618, which carry the digital data to and from 
computer 600, are exemplary forms of carrier waves transporting the information. 

[0119] Computer 600 can send messages and receive data, including program 
code, through the network(s), network link 620 and communication interface 618. In 
the Internet example, a server 630 might transmit a requested code for an application 
program tiirough Internet 628, ISP 626, local network 622 and communication 
interface 618. 
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[0120] The received code may be executed by processor 604 as it is received, 
and/or stored iii storage device 610, or other non-volatile storage for later execution. 
In this manner, computer 600 may obtain application code in the form of a carrier 
wave. 

[0121] It will be readily seen by one of ordinary skill in the art that the present 
invention fulfills all of the objects set forth above. After reading the foregoing 
specification, one of ordinary skill will be able to affect various changes, substitutions 
of equivalents and various other aspects of the invention as broadly disclosed herein. 
It is therefore intended that the protection granted hereon be limited only by the 
definition contained in the appended claims and equivalents thereof 
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